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Long term temporal correlations frequently appear at many levels of neural activity. We 
show that when such correlations appear in isolated neurons, they indicate the existence 
of slow underlying processes and lead to explicit conditions on the dynamics of these 
processes. Moreover, although these slow processes can potentially store information for 
long times, we demonstrate that this does not imply that the neuron possesses a long 
memory of its input, even if these processes are bidirectionally coupled with neuronal 
response. We derive these results for a broad class of biophysical neuron models, and 
then fit a specific model to recent experiments. The model reproduces the experimental 
results, exhibiting long term (days-long) correlations due to the interaction between slow 
variables and internal fluctuations. However, its memory of the input decays on a timescale 
of minutes. We suggest experiments to test these predictions directly. 

Keywords: neurons, temporal correlations, long memory, noise, input-output analysis, linear filters, power 
spectral density 



1. INTRODUCTION 

Long term temporal correlations, or "f^" statistics" (Keshner, 
1982), are ubiquitously found at multiple levels of brain and 
behavior (Ward and Greenwood, 2007, and refrences therein). For 
example,/^" statistics appear in human cognition (Gilden et al., 
1995; Repp, 2005), brain and network activity (measured using 
electroencephalograph or local field potentials Bedard et al., 2006, 
and refrences therin), and even Action Potentials (APs) generated 
by single neurons (Musha and Yamamoto, 1997; Gal et al., 2010). 
The presence of these long correlations in a neuron's AP responses 
suggests it is affected by processes with slow dynamics, which 
can retain information for long times. As a result, if these slow 
processes are also affected by APs, then the generation of each 
AP (indirectly) depends on a rather long history of the neuron's 
previous inputs and APs. This potentially allows a single neu- 
ron to perform complex computations over very long timescales. 
However, it remains unclear whether this type of computation 
indeed occurs. 

Cortical neurons indeed contain processes taking place 
on multiple timescales. Many types of ion channels are 
known, with a large range of kinetic rates (Channelpedia, 
http://channelpedia.epfl.ch/). Additional new sub-cellular kinetic 
processes are being discovered at an explosive rate (Bean, 2007; 
Sjostrom et al, 2008; Debanne et al., 2011). This variety is par- 
ticularly large for very slow processes (Marom, 2010). Such rich 
biophysical machinery can potentially modulate the generation of 
APs on long timescales. Evidence for such abilities was observed 
in recent works, which investigated how cortical neurons tempo- 
rally integrate noisy current stimuli (Lundstrom et al., 2008, 2010; 
Pozzorini et al., 2013). The temporal integration of the input was 
approximated using filters with power law decay, reflecting "long 
memory." However, these filters were fitted only up to a timescale 
of about 10 s (or equivalently, frequencies smaller than 10^' Hz), 



possibly due to the limited duration of the experiments, which 
involve intracellular recording. 

This raises the question - would the neuron still have long 
memory on timescales longer than 10 s? Generally, the answer 
may depend on the type of stimulus used. For example, certain 
ion channels may "remember" non-sparse inputs longer than 
sparse inputs (Soudry and Meir, 2010). Here, we focus on the 
case of the sparse (AP-like) input (Figure 1), imitating the "natu- 
ral" input for an axonal compartment which receives APs from a 
previous compartment. Such stimulation is used in various exper- 
iments (e.g., Grossman et al., 1979; De Col et al., 2008; Gal et al., 
2010). 

We find general conditions under which a neuron can gener- 
ate/^" statistics in its spiking activity, and show that this does not 
imply that a neuron has long memory of its history. Specifically, in 
order to generate/^" statistics slow processes should span a wide 
range of timescales with slower processes having a higher level of 
internally generated fluctuations (e.g., more "noisy," due to lower 
ion channel numbers). However, in a minimal model that gen- 
erates this behavior, slow processes do not retain memory of the 
input fluctuations beyond a finite "short" timescale, even though 
they are affected by the membrane's voltage. A main reason for 
this is that the "fastest adaptation process" in the model adjusts 
the neuron's response in such a way that any perturbation in the 
response is canceled out, before slower processes are affected. 

We fit the minimal model to the days-long experiments in 
Gal et al. (2010), where synaptically isolated individual neurons, 
from a rat cortical culture, were stimulated with extra-cellular 
sparse current pulses for an unprecedented duration of days. The 
neurons exhibited statistics, responding in a complex and 
irregular manner from seconds to days. The synaptic isolation 
of the neurons in the network, and their low cross-correlations 
indicate that these f^" fluctuations are internally generated in 
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FIGURE 1 I Stimulation Regime. (A) Stimulation consists of (extracellular) 
sparse current spikes, with inter-stimulus intervals Tm and Action Potential 
(AP) occurrences Vm- (B) An AP "occurred" if the voltage V crossed a 
threshold Vth following the (sparse) stimulus, with Tm » w. 



the neurons (Appendix D). We are able to reproduce their results 
(Figure 3), and predict that the neuron should remember pertur- 
bations in its input for about 10^ s (Figure 4). We suggest further 
experiments to test these predictions (Figure 5). 

The remainder of the paper is organized as follows. We begin 
in section 2.1 by presenting the basic setup. Then, in section 2.2, 
we present the general framework for biophysical modeling of 
neurons. Working in this framework, in section 2.3 we recall the 
mathematical formalism from Soudry and Meir (2014) and derive 
the power spectrum density for periodic input stimuli. Following 
a description of/^" behavior in section 3.1, we provide in section 
3.2 both general and "minimal" conditions for a neuron to display 
such scaling. In section 3.3 we consider the implications of the 
model for the input-output relation of the neuron, given general 
stationary inputs. In section 3.4 we demonstrate this numerically 
in a specific biophysical model which is fitted to the experimen- 
tal results of Gal et al. (2010). We conclude in section 4 with a 
summary and discussion of our results. An extensive appendix 
contains many of the technical details used throughout the paper 
(see Supplementary Material). 

2. METHODS 
2.1. PRELIMINARIES 

In our notation (•) is an ensemble average, i = a non- 

capital boldfaced letter x = (xi, . . . , x„)^ is a column vector 
[where (■)^ denotes transpose], and a boldfaced capital letter X 
is a matrix (with components X,„„). 

2.1.1. Stimulation 

As in Soudry and Meir (2014) we examine a single, synapti- 
cally isolated, excitable neuron under "spike" stimulation. In this 
stimulation regime, the stimulation current, /(f), consists of a 
train of identical short pulses arriving at times and ampli- 
tude /q. The intervals between the stimulation times are denoted 
Tm — tm+i — tm (Figure lA, top). We assume that the stimu- 
lation is sparse, i.e., r„, ^ tap, with tap being the timescale 
of an AP (Figure IB). Since the neuron is "excitable" it does 
not generate APs unless stimulated, as in Gal et al. (2010) 
(i.e., the neuron is neither oscillatory nor spontaneously fir- 
ing). However, after a stimulation the neuron can either respond 
with a detectable AP or not respond. We denote AP occurrences 
a Ym, where y,„ = 1 if an AP occurred immediately after the 
m-th stimulation, and 0 otherwise (Figure lA, bottom). Note also 
that Ym is not generally the same as the common count pro- 
cess generated from the APs by binning them into equally sized 



bins (Appendix B.l) - unless Tm is constant and equal to the 
bin size. 

2.1.2. Statistics 

We assume both Ym and Tm are wide-sense stationary (Papoulis 
andPillai, 1965). We denote = (y,„> to be the mean probability 
to generate an AP and T^, = ( r„,> as the mean stimulation period. 

Furthermore, we denote Ym = Ym — p* and f„, = Tm — T^, as 
the perturbations of y,„ and r„, from their means. An impor- 
tant tool in quantifying the statistics of signals is the power 
spectral density (PSD), namely the Fourier transform of the auto- 
covariance (Papoulis and Pillai, 1965). For analytical convenience, 
in this work we will use a PSD of the form 

oo 

SY{f) = T, J2 {YmYm + k)e-'^f^*'', (1) 

k=— oo 

with 0 r^' in Hertz frequency units. Note that this 

PSD is proportional to the PSD of the common binned AP 
(Equation 70), under periodical stimulus and for low frequen- 
cies - which is the regime under which we will investigate the PSD 
(similarly to the experiment Gal et al., 2010). We similarly define 
the PSD 

OO 

ST{f) = n J2 {fmfm + k)e-^^f^*'', (2) 

k=— oo 

and the cross-PSD 

oo 

Syr(/) = r, J2 {Ymfm+k)e-'-f^*'K (3) 

k=— oo 

2.2. GENERAL FRAMEWORK 

We model the neuron in the standard framework of biophys- 
ical neural models - i.e.. Conductance Based Models (CBMs). 
However, rather than focusing only on a specific model, we 
establish general results about a broad class of models. In this 
framework, the voltage dynamics of an isopotential neuron are 
determined by ion channels, protein pores which change confor- 
mations stochastically with voltage-dependent rates (Hille, 2001). 
On the population level, such dynamics are generically very well 
described by models of the form of Soudry and Meir (2014), 
(Equations 4-6) 

V=f{V,r,s,I(t)) (4) 
i = Ar(V)r-hr{V) + Br{V,r)^, (5) 
s = AAV)s-hs(V) + B,(V,s)^, (6) 

with voltage V, stimulation current /(f), rapid variables r 
(e.g., m, n, h in the Hodgkin-Huxley (HH) model Hodgkin and 
Huxley, 1952), slow "excitability" variables s e [0, 1]^ (e.g., slow 
sodium inactivation Chandler and Meves, 1970), white noise pro- 
cesses Ir^j (with zero mean and unit variance). Also, the matrices 
A,-/s and the vectors b^/s can be written explicitly using the kinetic 
rates of the ion channels, while the matrices Br/s can be written 
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using those rates in addition to ion channel numbers. Lastly, we 
denote 

as the diffusion matrices (Orio and Soudry, 2012). In these 
models the voltage and the rapid variables constitute the AP gen- 
eration, while the slow variables modulate the excitability of the 
cell. For simplicity, we assumed that r and s are not coupled 
directly, but this is non-essential (Soudry and Meir, 2014). The 
parameter space can be constrained (Soudry and Meir, 2012), 
since we consider here only excitable, non-oscillatory neurons 
which do not fire spontaneously^ and which have a single rest- 
ing state - as is common for isolated cortical cells, e.g.. Gal et al. 
(2010). 

2.3. THE POWER SPECTRAL DENSITY OF THE RESPONSE 

PSD-based estimators are central tools in quantifying long term 
correlations (Robinson, 2003; Lowen and Teich, 2005), and are 
commonly used in experimental settings - as in Gal et al. (2010). 
Therefore, in this section we focus on the PSD of the neural 
response under sparse stimulation regime (section 2.1) of a CBM 
(section 2.2). 

2.3. 1. Recap - previous mathematical results 

Typically, CBMs (Equations 4-6) contain many unknown param- 
eters, and are highly non-linear. Therefore, it is quite hard to fit 
them using a purely simulation based approach, especially over 
long timescales, where simulations are long and models have 
more unknown parameters. Therefore, we developed a reduc- 
tion method that simplifies analysis and enables fitting of such 
models. We refer the reader to Soudry and Meir (2014) for full 
mathematical details. 

In this method, we semi-analytically^ reduce the full model 
(Equations 4-6) to a simplified model, under the assumption 
that the timescales of rapid and slow variables are well separated. 
Given another assumption, that the neuron dynamics are suf- 
ficiently "noisy," we can linearize the model dynamics, so that 
(Soudry and Meir, 2014, Equation 12) 

Ym = (s (fm) - s*) -I- e™, (7) 

where is a white-noise signal with zero mean and variance 

= p^, — (recall p^, is the mean probability to generate an 
AP) and s^, (the excitability fixed point) and wj (an "effective 
weight" of component Sj) can be found self-consistently together 
withp^, as a function of (SoudryandMeir, 2014, Equation 10). 
After these quantities are found, an expression for the output PSD 
Sy if) in this model can be written explicitly. We let X+,X- and 
Xq denote the averages of the quantity X^ during an AP response, 
a failed AP response and rest, respectively. Also, we denote 

X, = TAPT-I {p,X+ + (1 - + (1 - TApT-')Xo 



I.e., if Vt:7(t) = 0, then the probabihty that a neuron wiU fire is 
negligible - on any relevant finite time interval (e.g., minutes or days). 

semi-analytic derivation is an analytic derivation in which some terms are 
obtained by relatively simple numerics. 



as the steady state mean value of Xj [this would be X (p^, , T^,) in 
Equation 7 in Soudry and Meir, 2014]. For example, A^, and D^, 
are the respective steady state means of Aj and Dj. Additionally, 
we denote (definition below Equation 12 in Soudry and Meir, 
2014) 

a = rAp((A+ - A_)s, - (b+ - b_)) (8) 

as a "feedback" vector (see Figure IC in Soudry and Meir (2014) 
to understand this interpretation), and (Soudry and Meir, 2014, 
Equation 14) 

H, (f) = (27r^I - A, - T-'aw^y (9) 

as the "closed loop transfer function" (including the effect of 
the feedback), with I being the identity matrix. Using the above 
notation, we can derive the PSD of the response. Given a peri- 
odical stimulation (T^ = 0) we obtain (Soudry and Meir, 2014, 
Equation 13) 

Sr(f) =w^H, (-/)D,H7(f)w 

+ r,ff2 1 i+r-iw^H, (-/) . (10) 

Though Equation (10) relies on two simplifying assump- 
tions, extensive numerical simulations (Soudry and Meir, 2014, 
Figures 3-5) showed that this expression is rather robust and 
remains accurate in many cases even if these assumptions do not 
hold. Therefore, in this work we will always assume that Equation 
(10) is accurate. 

2.3.2. The effect of feedhack 

In the neuron, the slow excitability variables s affect the response 
of the neuron, which, in turn, affects the dynamics of the slow 
excitability variables. To simplify analysis, it is desirable to "iso- 
late" this feedback effect. In order to do this, we apply the 
Sherman Morrison lemma to Equation (9), 

w^H, (-/) = wTh„ (-/) (l - r-iwTH„ (-/) , 
with 

H„{f) = {27tfil-A,)-' (11) 

being the "open loop" version of (f) (i.e., if a is set to zero). 
Using this in Equation (10) we obtain 

SY(f)=S"y(f)\K(f)\-\ (12) 

with Sy (f) being the "open loop" version of Sy (/) (i.e., Sy (/) 
with a set to zero), 

Sy if) = + wTh„ (-/) D,hJ if) w (13) 

and K (/) determines the effect of the feedback 

^(f) = l-r-'wTH„ (-/)«. (14) 
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Note that k (f) depends on the feedback through the variable 
a. If a— 0, for example, the kinetic rates of s are not sensi- 
tive to AP occurrences ^. In that case {/) — ^ 1 and Sy (f) — > 

Syif)- 

2.3.3. Partial fractions decomposition 

In order to simplif)f analysis, we decompose the vector expressions 
in Equations (13, 14) to partial fractions. 

If A^, is diagonalizable, than we can write Equation (13) as 
(Appendix A. 1 ) 



M 



[{2nff + Xl' 



(15) 



where the poles Xjt are the inverse timescales of the slow variables 
(the eigenvalues of A,, ) , arranged from large to small according to 
their magnitudes (0 < \Xm\ < Um-iI < ••• < |Ai|)and 



M 



Ck 



^ WkDkjWj 



2Xk 



(16) 



being the amplitude of these poles, with Dj^j and Wk being the 
respective components of D^, and w in a basis in which A^,, is 
diagonal. Note that, Vfc, Re [kk] < 0 (from the properties of A^,). 
Using a similar derivation for k (f), we obtain 



M 



-(f) = i-E 



i:=l 



WkUk 

Infi — kj^ 



(17) 



with Uk and Wk being the respective components of a and w in a 
base in which A^ is diagonal. 

2.3.4. Example -a "diagonal" model 

For concreteness, we demonstrate our results on a simple model 
in which A^, is a diagonal matrix and, as a result, (which 
depends on A* ) is also diagonal. In this "diagonal" model all the 
components of s are uncoupled (i.e., belong to different channel 
types). Equation (6) can be written as (Soudry and Meir, 2014, 
section 4.1) 

h = Sk (V) (1 - Sk) - Yk (V) Sk + (Ts,k (V, Sk) ?a (18) 



Vfc e {1, . . . , M}, where aa (V, s)=[(Sk (V) (1 - Sk) + yk (V) Sk) 

_,-|l/2 

N J. and are the number of slow ion channels of type k. 

Similarly as before, K+.j:, y-,k and yo,J: denote the averages of the 
kinetic rate yk (V) during an AP response, a failed AP response 
and rest, respectively. In addition 

y*,k = TApT^' {p*y+,k + (i -p*) y-.k) + (i - tapT^^) yo,k 



For example, this can happen if the kinetic rates aU have low voltage 
threshold, resulting in A-|- ~ A_ and b-|- ~ b_ . 



is the average yk (V) in steady state. We use a similar notation for 
S. Therefore 

(■^*)kk = ~y*.k — S^:,k 

1 y*,kS*,k 



kk 



Ns.k y*,k + i5*,i 

with zero on aU other (non-diagonal) components and 

{r*,k {s+.k - S-,k) - {y+.k - y-.k) s^,,k) 

ak = Tap — ; • 

y*.k + o*,k 

Therefore, in Equations (15, 16) we have. 



kk 



Ck 



-y*,k — Si,,k 



wlDkk ■■ 



^k Y*.k&*.k 
Ns.k y*,k + 



(19) 

(20) 
(21) 



Importantly, by tuning the parameters M, yk (V) , Sk (V), Njjt 
and Wk we seem to have complete freedom in determining Aj-, 
Ck and ak (Equations 19-21). This, in turn, would give complete 
freedom in tuning Sy (f) and k (f). Therefore, it seems that for 
any CBM (i.e., not only diagonal models) we can find an equiv- 
alent diagonal model - which produces exactly the same PSD of 
the response. 

The only caveat in the previous argument is that in non- 
diagonal models kk can be complex, but not in a diagonal model, 
since the kinetic rates yj^ (V) and 5jt (V) must be real numbers. 
How would the situation change if some of the poles had com- 
plex values? Complex poles (i.e., for which Im [kk] > 0) always 
come in conjugate pairs. These pairs behave asymptotically (i.e., 
for Inf ^ or 2nf <^ \kk\) very similarly to two real poles, 
with an additional "resonance" (either a bump or depression) in 
a narrow range in the vicinity of these poles (i.e., 27r/ ~ \kk\) (see 
Appendix A.2, or Oppenheim et al, 1983). 

3. RESULTS 

3.1. BACKGROUND ON f'" STATISTICS 

As observed in Gal et al. (2010), the responses of isolated neu- 
rons exhibit long-term correlations robustly under sparse pulse 
stimulation (Figure 1 and section 2.1). Signals with such long- 
term correlations are often described by the term "f^" noise." 
This is because the Power Spectral Density (PSD, Papoulis and 
Pillai, 1965) is a central tool in detecting and quantifying such 
signals (Robinson, 2003; Lowen and Teich, 2005). As the name 
implies, if the AP pattern is a "f^" noise signal" then its PSD 
(Equation 10) has a/^" shape 



SY{f)<xf-\ 



(22) 



where the PSD is defined here as in Equation (1). As is usually 
the case for most f^" phenomena, Equation (22) is true only 
in a certain range fmm Sf< /max, and with 0 < a < 2. Note 

''i.e., in all neurons for which 0 < < 1. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



April 2014 I Volume 8 1 Article 35 | 4 



Soudry and Meir 



Long correlations without long memory 



also that if a > 1, then /min > 0 necessarily^. Such f^" behav- 
ior is considered interesting due to its "scale-free" properties, 
which can sometimes indicate a "long memory," as explained in 
the introduction. Therefore, it is interesting to ask the following 
questions: 

1. What is the biophysical origin of the/^" behavior? 

2. Does this/^" behavior imply that the neuron "remembers" its 
history on very long timescales (hours and days)? 

We aim to answer the first question in section 3.2, focusing on the 
case of periodical stimulation Tm = T^,, as in Gal et al. (2010). 
The second question is addressed in section 3.3, where we exam- 
ine a general sparse stimulation process T„,. Finally, in section 
3.4.2 we fit a specific CBM (which is an extension of a previous 
CBM) so it adheres to this set of minimal constraints. We numer- 
ically reproduce the experimental results of Gal et al. (2010) and 
demonstrate our predictions. 

3.2. BIOPHYSICAL MODELING OF f" STATISTICS 

As we explained in the introduction, neurons contain a large vari- 
ety of processes operating on slow timescales. These processes 
are, in many cases, not well characterized or contain unknown 
parameters. Therefore, it is hard to model the behavior of the neu- 
ron on slow timescales with a CBM using only simulation. With 
so many unknowns, an exhaustive parameter search is unfeasi- 
ble ^. Fortunately, since we derived a semi-analytic expression for 
the PSD (Equation 10), starting from some initial "guess" (as to 
which process to include, and with what parameters), it is rel- 
atively straightforward to tune the parameters so that the CBM 
reproduces the experimental results (i.e., by maximizing some 
"goodness of fit" measure). 

However, even if a specific model could be found to repro- 
duce the experimental results, it would stiU be unclear whether 
or not this is would be a "useful" model - one which can be used 
to infer the biophysical properties of the neuron, or its response 
to untested inputs. The first problem is that CBMs are highly 
degenerate, where different parameter values can generate simi- 
lar behaviors^, so we can never be sure if the "correct" model was 
inferred. The second problem is that it is unclear whether a "cor- 
rect" model would be generally useful - since different neurons 
from the same type can have very different parameters (Marder 
and Goaillard, 2006). 

In order to address the first problem, initially, in section 3.2.1 
we analyze Equation (10), and attempt to answer a more gen- 
eral question - what class of CBM models can generate the 
experimental results? We find "rather general" sufficient con- 
ditions - i.e., which, given a few assumptions, also become 
necessary conditions. Next, in section 3.2.2, we aim to find a 
"minimal" set of constraint on a CBM to fulfill theses conditions. 



Qualitatively, these conditions indicate that, in order to reproduce 
the experimental results, a general CBM must: 

1. Include only a finite number of ion channels of each type 
(implying a stochastic model). 

2. Include few slow processes with timescales "covering" the 
range of timescales over which Sy (/) oc is observed. 

3. Obey a certain scaling relation (with an exponent of 1 — a), 
implying that slower processes are more "noisy." 

More detailed explanations of these conditions, and a concrete 
example, are provided in the following two subsections. 

3.2. 1. General conditions for f~" statistics 

In this section we derive general conditions on the parameters 
of a CBM (section 2.2) so it can generate robust f^" statistics 
in Sy (f). Here, we focus on the case of sparse periodical input 
Tm = T^: ^ Tap (as in Gal et al, 2010). 

This analysis is based on the decomposition of the PSD Sy (/) 

as a ratio of Sy (/) and the feedback term | k (f)\ . Recall that 
Sy (f) OC f^" is robustly observed for all stimulation parame- 
ters - even when p^, is near 0 or 1 (see section 3.1). Note that 
one can arbitrarily vary p^, by changing the stimulation param- 
eters (such as 7o or T*). It is straightforward to show that when 
p* ^ 0 or p^, ^ 1, the effect of feedback is negligible ^ and 
therefore Sy (f) ~ Sy (/). This implies that, at least for some 
simulation parameters, Sy (f) cx/^". For this reason, and for 
the sake of analytical simplicity, we first develop general condi- 
tions so that Sy (f) oc f^", and later we discuss the effects of the 
feedback a: (/). 

Note from Equation (15) that if M (the dimension of s - 
the number of slow processes) is finite, one can have Sy (/) cx 
f^" exactly if and only if a = 0 or 2. However, these values 
are far from what was measured experimentally (Equation 42). 
Therefore, Sy (f) cx can be generated exactly only in some 
limit (in which M is infinite), or approximately (if M is 
finite). Also, note that if 27r/ » then S"y (f) - T^aj oc 
/^^. Additionally, if 27r/ ^ |XmI> we have Sy (/) ~ constant. 
Therefore, Equation (15) can generate Sy (/) oc/^" with 0 < 
a < 2 only for \Xm\ < 2jr/ < \Xi\. 

Next, we explain when this becomes possible. For simplicity 
assume that in Equation (15) T^^cTj? is negligible and all the poles 
are real (the effect of complex poles will be discussed below). We 
define the following pole density 



M 



(23) 



where S (•) is Dirac's delta function. Using Equations (15 and 23) 
we obtain 



^Otherwise, 0.25 > p« - = |y^| > 2 /^"" Sy (f) (i/ = oo, which is a 
contradiction. 

^Also simulations take a long time, since experiments, as in Gal et al. (2010), 
are days-long. 

^E.g., in Equation (16) many different parameters would give the same cj. 



S^yif)=f 



(24) 



Near the edges, w - 
.(f) 



0 (Equation 80 in Soudry and Meir, 2014), and so 
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For I^mI <^ Inf <^ \ki | and 0 < a < 2, Equation (24) becomes 
Sy if) = Cf-" (25) 
if and only if (Appendix A. 3) 

p(k) = po\X\'-" (26) 

in the range \Xi \ > \X\ > |>.mI, with po = 27r^^Csin (;ra/2). 
Therefore, p (X), the distribution of the poles, must scale similarly 
to Sy (f) (but with a different exponent). 
Several comments are in order at this point. 

1. It was previously known that, in a linear system, a f^" 
PSD could be generated using a similarly scaled sum of real 
poles (Keshner, 1979, 1982). The novelty here is two-fold: 
(1) Quantitatively analyzing the PSD of CBMs (which are 
highly non-linear) in a similar way (through Equation 10) (2) 
Finding that condition 26 is not only sufficient, but necessary. 

2. Formally, Equation (26) can be exact only in the continuum 
limit where the number poles is infinite and they are closely 
packed. However, in practice. Equation (25) remains a rather 
accurate approximation even if the poles are few and well sep- 
arated (Figure 2A), as we shall demonstrate in the next section 
(as in Keshner, 1979, 1982). Clearly, for simulation purposes, it 
is beneficial to use a CBM with a finite number of (preferably, 
few) poles . 

3. We have assumed that all the poles are real. What happens if 
some of the poles are complex? Recall (section 2.3.4) that if 
some poles have complex values then Sy (/) also has "reso- 
nances" (bumps or depressions) in a narrow range near these 
poles. Technically, scaling these resonance peaks can also be 
used to approximate Equation (25) (Figure 2B). However, we 
did not pursue this method here since it would require signif- 
icantly more poles and would be much harder to implement. 

4. Note that so far we have discussed only Sy (/). One can per- 
form a similar analysis directly onSy (f). However, we find it is 
easier to first simpliiy k (/) and then use Equation (12). From 



Equation (12) the PSD Sy (f) will have a power-law shape in 
the range {XmI <SC 27r/ <SC I'^il if, in that range: either (1) the 
magnitude of k (f) is constant or slowly varying, or (2) k (f) 
also has a power-law shape. In the first case the exponent of 
Sy (f) will be the same as the exponent of Sy (f), and in the 
second case the exponent will differ. The conditions for both 
cases can be derived similarly to our analysis of Sy {/) . We 
demonstrate this next, in a more specific context. 

3.2.2. A minimal model for f~" statistics 

In the previous section we found general conditions under which 
Equation (13) gives Sy (f) oc/^". In this section, we aim is to 
generate Sy (/) oc over/min < / < /max in a minimal model, 
in which M (the dimension of s) is as small as possible. As 
explained in section 2.3.4 we do not lose any relevant gener- 
ality if we restrict ourselves to the case where A,, is diagonal 
(Equation 18). From Equation (26), we know that \kk\ must 
"cover" the frequency range /min < / < /max- In order for M to be 
small, we choose kk to be uniform over a logarithmic scale (sim- 
ilarly to Keshner, 1982), so Xk oc 6*^ with 6 < 1. The "simplest" 
way to achieve this is to have (see Equation 18) 

yk(V) = n{V)€^-' ; SkiV) = SiiV)e^-' (27) 

so 

kk = kie''-\ (28) 

In order for ki^/ (Itt) to cover the range [/min,/max] we require 
that 

l^ll > 2;r/max; I^mI = l^ll 6^"^ « 27r/min. (29) 

Given M, this sets a constraint on 6. In order to have scaling in 
p (k), as in Equation (26), we also require that cj^ oc l^l^^" dk oc 
^(2-a)k^ since dk = k^ — kk- i oc e''. Therefore, from Equations 
(21) and (20) we have 

so that Sy (/) oc f^" . Therefore, we require that Wk oc e^'^'^, 
Ns,k oc 6"''' with Ifi + V = a — i. For /x > 0 the slower processes 
(larger k) have larger weight. For v > 0 slower processes have 
a smaller number of ion channels (therefore, they are more 
"noisy"). 

In Appendix A.4, we investigate what type of scaling will 
generate also Sy (/) oc f^" , taking into account the effects of 
feedback (through k (/)). We conclude that, because of the feed- 
back, a value of > 0 would not change the exponent of Sy (f) 
over a "reasonable" range of parameters (i.e., assuming v > —2). 
Therefore, the simplest way to generate Sy (/) oc f^" would be to 
take /X = 0. In this case, we have (Equation 59), for — 1 < v < 1 
andl^Ml «27r/ «; 

1 f-(l+v) 

Sy (/)oc ■' — = — , (30) 



log «?•(/) 




FIGURE 2 I Generating f" PSD using a finite number of poles - a 
graphic description. Using partial fraction decomposition (Equation 15) 
S° (f) (X f^" (blue) can be approximated (on a log-log scale) in two distinct 
ways: (A) Using a sum of a real poles (green), with scaled amplitudes 
(approximating Equation 26) (B) Using a sum of complex poles (orange), 
with scaled "resonance peaks" (Equation 46). In this work we focus on the 
first case (A), since it is simpler and requires far fewer poles. 
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FIGURE 3 I The measures of "scale free" rate dynamics in the 
HHMS model - comparison of the experimental data from Gal 
et al. (2010) and a simulation of the extended HHS model (solid 
and dashed lines, respectively). We use here the same measures as 
in Figure 6 in Gal et al. (2010): (A) The firing rate fluctuations 
estimated using bins of different sizes {T = W, 30, 100, and 300s) and 
plotted on a normalized time axis (units in number of bins), after 



subtracting the mean of each series. (B) CV of the bin counts, as a 
function of bin size, plotted on a log-linear axis. (C) Firing rate 
periodogram. (D) Detrended fluctuations analysis. (E) Fano factor (FF) 
curve. (F) Allan factor (AF) curve. (G) Length distribution of 
spike-response sequences, on a half-logarithmic axes. (H) Length 
distribution of no-spike-response sequences, on a double-logarithmic 
axes. For additional details on measures used see Appendix B.I. 



where the logarithmic correction arises from the effect of feed- 
back /f (/). A few comments on Equation (30) are in order at this 
point. 

1. Due to the logarithmic correction, in order to approximate 
Sy (f) oc f^" it is a reasonable choice to set v slightly higher 
than a — 1, e.g., 

v = a-Q.9. (31) 

2. Even if there is no scaling in the parameters (i.e., /x = v = 0), 
we obtain Sy (f) oc/^' (neglecting logarithmic factors). 

3. Equation (30) is based on asymptotic derivation, which is 
correct in two opposing limits ("sparse" or "dense" poles, 
Appendbc A.5), indicating that these results are rather robust 
to parameter perturbations. 

4. The magnitude of the ion channels number Nj i is inversely 
proportional to the magnitude of Sy (/) (i.e., its proportion- 
ality constant), while the value Wi (the magnitude of the 
weights) does not affect Sy (f). 

5. When Njj ^ oo we have Sy (/) 0, implying that in the 
deterministic limit, such a CBM does not generate f^" noise 
(in accordance with our results from Soudry and Meir, 2012). 



3.3. THE INPUT-OUTPUT RELATION OF THE NEURON 

In the previous section we derived minimal biophysical con- 
straints under which a neuron may generate f^" statistics in 
response to periodic stimulation. In this section we explore the 
input-output relation of the neuron under these constraints, in 
the case where the inter- stimulus intervals T„, form a general 
(sparse) random process. We decompose the neuronal response 
into contributions from its "long" history of internal fluctuations 
and its "short" history of inputs, quantifying neuronal memory. 

3.3. 1. The linearized input-output relation 

Recall that f„, = T,„ - T^, with T^, = {!„,) and Sj if) the PSD 
of Tm (Equation 2). As explained in Soudry and Meir (2014), for 
a general CBM' we can decompose Y,„, the fluctuations in the 
neuronal response, to a linear sum of the history of the input and 
internal noise, i.e., 

oo oo 
i"m = ^ hf^^Tm-k + ^ hf^Zm-k, (32) 
i: = 0 k = 0 



'i.e., Equations (4-6), with the same assumptions as we had in section 2.3.1. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



April 2014 I Volume 8 | Article 35 | 7 



Soudry and Meir 



Long correlations without long memory 




FIGURE 4 I System decomposition into external (input - r(f)) and intemal (fluctuation - z (f)) filters. For a fitted HHMS model, H""' (f) is a low pass 
filter with cutoff < 10"^ Hz while Hint (0 ~ f'"'^ for f < 10"^ Hz. 
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FIGURE 5 I Input memory in fitted model. (A) Comparison of |Sy7 (f)| of 
the fitted model ("Model") to that estimated from the experimental 
confirms ("Experiment") the prediction of the input filter IH""' for 
probed range. (B) This filter ("Approx") can be probed more accurately by 
peaks of Y{f) ("Simulation"), by applying a "sum of sines" input 
(Equation 43). 



with the filter /j™' used to integrate external fluctuations in the 
inputs, and the filter used to integrate z„,, a zero mean and 
unit variance white noise representing internal fluctuations (e.g., 
ion channel noise). It is easier to analyze this I/O in the fi'equency 
domain, where Equation (32) becomes (Soudry and Meir, 2014, 
Equation 20) 



y(/)=H«t(/)r(f)+H-'{/)z(f). 



(33) 



where we define to be the Fourier transform of X{t). 

Together, H'^''* (/) and H'"' (f) describe the Y,„ neuronal 

I/O at very long timescales. 

Note that these filters are related to the PSDs, in the following 
way 

Syt if) = H'"^ (-/) St if) , (34) 
Sy if) = (/) I' St (f) + |h'"' (/) I' (35) 



3.3.2. The shape of the input-output filters 

For a general CBM, we can derive semi-analytically the exact form 
of the filters in Equation (33) from its parameters, as we did 
for Sy (f). For example, if = 0 (periodical input), then also 
St (f) = 0, and so 



\Hint(f)\=SY(f). 



(36) 



where Sy (/) is the PSD we derived previously (Equation 10). 
Additionally, we obtain (Soudry and Meir, 2014, Equation 17) 



with 



H^'"(f) = r-iw'H,(f)d 



d = Aqs^, — bo. 



(37) 



(38) 



Next, we find both filters for the minimal model described in 
section 3.2.2. Recall that in this model 



k J k 



(39) 



with ai and di, respectively given by Equations (8) and (38). To 
simplify analysis, we derive an asymptotic form for both filters, 
for the cases \Xm\ ^ 2nf <JC and lirf ^ \Xi \ ■ First, from 
Equation (36), and Equation (59), we find 



Hi, 



jr"/Vln/ ,if |XMl«27r/«|Ai| ^^^^ 
' [constant ,if27r/»|Xi| 



Similarly, from Equation (37), we find (Appendix A.6) that for 
the minimal model the interpolation between the two asymptotic 
cases is monotonic, so we can approximate 



qdi 

Ijrfi — qai 



(41) 



where we recall that Syr if) is the cross-PSD (Equation 3). 
Notably, from Equation (35), if the input to the neuron is 
not periodical (so, St (/) / 0), then the PSD Sy (/) should be 
the same as calculated previously, except for the addition of 

|H-*(f)|'Sr(/). 



where q = (I — e) ^ ^wi. A few comments on Equations (40, 
41) are in order at this point. 



1. We found that H™' 
qai/2jT while H'"* 



(/) is a low pass filter with a pole at/ext = 
(y) ~ f-oin for 27r/ «; |Ai |. Consequently 
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in the temporal domain (Equation 32), for large t (i.e., large k), 
the neuron's memory of its external input decays exponentially 
(/jcxt ^ e^f^'^*^), while its memory of its internal fluctuations 
decays as a power law (/zj,"' ~ kr''^^"!^'*). Therefore, the input 
memory has a finite timescale (equal to/g^/), while the mem- 
ory of internal fluctuations is "long" (with a cutoff only near 

2. It is perhaps surprising that Equation (37), which has mul- 
tiple poles, becomes a low pass filter with a single pole /i . 
The derivation (Appendix A.6) gives two main reasons for 
this. First, the scaling of Wk and in Equation (39) induces 
only a weak (logarithmic) scaling of the poles in open-loop. 
Second, even this weak scaling is canceled by the effects of the 
feedback. 

3. Naturally, other models may have a different shape of H™* (/) . 
This could be probed directly, as we explain later, in section 
3.4.3. 

3.4. MODELING EXPERIMENTAL RESULTS 

In this section we apply our results to experimental data, 
described in section 3.4.1. In section 3.4.2 we implement the set 
of "minimal constraints" we found in section 3.2.2 in a specific 
CBM, and fit it to experimental data in which Sy (f) oc f^" . The 
analytical results in section 3.2 suggest that this specific CBM 
is a "reasonable" representative of the family of CRMs that can 
generate the experimental results. Other members of this fam- 
ily can be reached by varying the parameters within the (either 
minimal or general) constraints. Next, in section 3.4.3 we use 
our results from section 3.3.2 on the fitted model. We show that, 
although internal fluctuations in the model can affect the neural 
response on a timescale of days, the memory of the input is only 
retained for a duration of minutes. We suggest specific experi- 
ments to test this prediction. In section 3.4.4 we suggest further 
predictions 

3.4. 1. Experimental details 

The experiment from Gal et al. (2010), where a single synapti- 
cally isolated neuron, residing in a culture of rat cortical neurons, 
is stimulated periodically with a train of extracellular short cur- 
rent pulses with constant amplitude 7o- The observed neuronal 
response was characterized by different modes (Gal et al., 2010, 
Figure 2). We focus on the "intermittent mode" steady state, in 
which 0 < p^, < 1 (i.e., sometimes the stimulation evokes an AP, 
and sometimes it does not). The patterns observed in Ym, the 
AP occurrences timeseries, are rather irregular (Gal et al., 2010, 
Figure 2E), span multiple timescales (Gal et al, 2010, Figures) 
and variable (i.e., patterns are not repeatable Gal et al., 2010, 
Figure 9A). More quantitatively, as indicated by the analysis (Gal 
et al., 2010, Figure 6), for all intermittently firing neurons, the pat- 
terns in Ym fall into the category of "/^" noise" where the value 
of a varied significantly between neurons - with 

a = 1.43 ± 0.35 (42) 

(mean ± SD). As we explained in section 3.1, this/^" behavior is 
true only in some limited range /min < f < /max- From the exper- 
imental data, (Figure 6C in Gal et al., 2010) it can be estimated 



that /min < 10 ^Hz and /max ~ 10 ^Hz. Also, since a > 1, then 
0 < /min (see section 3.1). 

3.4.2. The HHMS model - a biophysical implementation of the 
minimal constraints 

In our previous work (Soudry and Meir, 2012) we already fitted 
a model that fits many of the "mean" properties of the neu- 
ronal response (e.g., firing modes, transients and firing rate). This 
model is an extension of the original Hodgkin-Huxley model 
which includes Slow sodium inactivation (Chandler and Meves, 
1970; Fleidervish et al, 1996) (The HHS model, see Appendk 
C.l). In order to maintain this fit with the experimental results, 
we extend the HHS model with additional slow components, 
obeying Equation (18). We denote this as the HHMS model 
(Hodgkin Huxley model with Many Slow variables. Appendix 
C.2). The equations are identical to the HHS model, except 
that in the voltage equation (Equation 73) g]VaS is replaced by 
SNaM^^Ylh=i^i' where si has the same equation as s in the 
HHS model (Equation 77). By symmetry, this gives identical 
weights to component s,- (i.e., Vfc : wk = wi). The remaining 
rates (for k > 2) are chosen according to our constraints, so 
Yk (V) = y (V) 6*=- K 4 (V) = S (V) e''- ' (as in Equation 27), 
where y (Y) and 8 {V) are taken from the HHS model (Equation 
77) and also N^^k = N^e"'^. Therefore, the only free parameters are 
e, M, Ns, V and 7o (-fo is the current amplitude of the stimulation 
pulses). 

This model can be used to fit the experimental results for 
any a e [0, 2). We performed a numerical simulation of the full 
equations (Equations 4-6) of the HHMS model under periodi- 
cal stimulation with T^, = 50 ms. We aimed to fit an experiment 
from Gal et al. (2010), which had a similar stimulation and 
exhibited Sy (f) ~/ "> with a = 1.4 (which is approximately 
the average a value measured in Gal et al., 2010). The cur- 
rent amplitude Iq was set to Iq = 7.7 IJ,A so that the model 
would have the same mean response probability p^, ~ 0.4 as in 
the experimental data (using the self consistent equations for 
pit from Soudry and Meir, 2014). We chose M = 5 and e = 
0.2 in order to satisfy constraint Equation (29) with a mini- 
mal M. We chose v = 0.5 to satisfy Equation (31). Lastly, we 
chose Ns = 10* in order to fit the magnitude of the Sy {/). 
This reproduced all the scaling relations observed experimentally 
(Figure 3). 

3.4.3. Predictions - probing the input-output relation of the neuron 

After fitting the HHMS model to the experimental results, we can 
examine its resulting linearized input-output relation, described 
by the filters H'""' (/) and H'"* (/) (Equation 33). The H'"' (/) 
filter integrates internal fluctuations , while the H'^^ (/) filter 
determines how external fluctuations (in the input) affect its 
response. 

In accordance with the asymptotic forms in Equations (40) 
and (41), we find that H™* (/) is a low pass filter with a pole/ext ~ 
10^2 Hz (Figure 4, green) while H'"' (f ) ~ /""/^ for /^in < 
/ < 10^^ Hz (Figure 4, red) with /min < 10^^ Hz. Therefore, as 
explained in section 3.3.2 this model implies that the response of 
the neuron is affected by internal fluctuations over the scale of 
days (~ /mill) or more, generating the f^" behavior we observe 
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in Figure 3. However, external input is remembered only for 
minutes (~/ex/). 

Next, we examine two methods which allow us to probe 
H™' (f ) directly and examine these predictions. 

First, a simple method to probe the external input filter 
H™' (/) is through Equation (34). Allowing reliable estimation 
of H™* (/) in a certain frequency range requires a random pro- 
cess stimulation for which |H™* {f)\^ St (/) » |H'"' (f) \ in that 
range, as explained in Appendix B.2. To demonstrate this method 
we estimate Sty (f) from the existing experimental data taken 
from Gal and Marom (2013), in which St (/) ~ /"'^ (above some 
lower cutoff). In Figure 5A we compare this estimation with 
Sty (f) in the HHMS model in a limited range where St (f) is suf- 
ficiently high for estimation to be accurate. It is similar to Sty (/) 
from our fitted HHMS model, validating our estimate of H™' (/) 
for that frequency range. 

Second, The filter H™' (/) could be probed more accurately 
and at lower frequencies - by sinusoidaUy modulating the input 
(the internal-stimulus intervals), analogously to the sinusoidaUy 
modulated input current used in Lundstrom et al. (2008, 2010) 
and Pozzorini et al. (2013), 

I 

Tm = Tamp ^ sin [infiT^^m) . (43) 
;= 1 

As we explain in Appendix B.3, in this case the output of the 
neuron would be 

L 

Ym = J2 l"'"" (/')hin (27r/,r,m + IH'"^ (f,)) + "noise." 
(= 1 

This allows us to easily estimate |H™*(/)| using the peaks of 
T^JipY (f) (the Fourier transform of T^J^pYm) at frequencies 
/;, as we demonstrate in Figure 5B, using our fitted HHMS 
model. 

3.4.4. Additional predictions 

As explained in Gal et al. (2010) and Soudry and Meir (2012) the 
latency of the AP can serve as an indicator of the cell's excitabil- 
ity. Specifically, this is true in the HHMS model, for periodical 
stimulus and = 1, where the PSD of the latency, Sl (f), is 
a shifted and scaled version of Sy (f) with ^ 1 (see section 
4.4.6 in Soudry and Meir, 2014). Therefore, in the HHMS model 
we also have Si (f) a f^" approximately (neglecting logarithmic 
factors). 

Next, suppose we vary some measurable stimulation parame- 
ter, such as the mean stimulation rate ' . How would this affect 
the shape of the filters we derived? The analytical results allow us 
to calculate this explicitly in the HHMS model. 

First, we consider the gain of the external input filter H'^^^ (/) 
(i.e., H'^^ (0)). As we explain in Appendix A.7, if/ <^ /cutoff' than 

H^'^{/)«p*r-l=/out, (44) 

which is the mean firing rate of the neuron - a simply measurable 
quantity. 



Second, how would Hint if) change if T^, is varied? Since 
Hint (/) is directly measurable only through Sy (/) (Equation 36), 
we consider Sy (f) instead. From Equation (59) it is clear that if 
Sy (f) ~ f^" approximately at low frequencies then the exponent 
a should not depend much on any external parameter (assuming 
0 < pit < 1). This was observed experimentally when the stimu- 
lation rate (T^^) was varied, as can be seen in Figure IG in Gal 
and Marom (2013). 

4. DISCUSSION 

4.1. GENERATING PSD 

In this work we aim to explain biophysicaUy the phenomenon of 
/^^ behavior in the response of isolated neurons, and explore 
its implications on the input-output relation of the neuron. We 
do this under a regime of sparse stimulation (Figure 1), and in 
the biophysical framework of stochastic conductance-based mod- 
els (CBMs, Equations 4-6). In this setting our analytical results 
(Soudry and Meir, 2014) can be used to derive a closed form 
expression for the Power Spectral Density (PSD, Equation 10) 
based on the parameters of the slow kinetics in the CBM. This 
PSD is affected by the closed-loop interaction - the slow dynamics 
affect the AP response, which, in turn, feeds back and affects the 
kinetics of the slow processes (section 2.3.2). Moreover, the con- 
tribution of each slow process to the PSD can be exactly quantified 
(section 2.3.3), as we demonstrate using a simple model (section 
2.3.4). 

These mathematical results expose the large parameter degen- 
eracy of CBMs (Marder and Goaillard, 2006; Soudry and Meir, 
2014), i.e., that many "different" models will quantitatively pro- 
duce the same behavior. Due to the degeneracy of CBMs, we first 
aimed to derive rather general sufficient conditions for the gen- 
eration of/^" noise in a CBM (section 3.2.1). These conditions 
indicate which types of CBMs can generate the observed behavior. 
We show that, in order to generate /^" behavior, neurons should 
have intrinsic fluctuations (e.g., due to ion channel noise), and 
have a number of slow processes with a large range of timescales, 
"covering" the entire range over which f^" statistics is observed. 
Furthermore, the parameters of these processes must be scaled 
in a certain way in order to generate f^" noise with a specific a 
(Equation 26). 

We implement these constraints in a minimal CBM (section 
3.2.2), in which the slow processes are uncoupled, except through 
the voltage, as in Soudry and Meir (2012). Initially, we find that 
the specific scaling relation can be achieved either by scaling the 
(1) conductances or (2) the ion channel numbers. This scaling 
implies that slower processes will be either (1) "stronger" or (2) 
"noisier." However, the "feedback" effect in the model (the slow 
process being affected by the APs) prevents f^" statistics from 
being generated in case (1). In contrast, option (2) can robustly 
generate the observed f^" statistics in the neuronal response for 
any 0 < a < 2 (Equation 30 and Figure 6). 

Naturally, outside of the framework of CBMs (Equations 
4-6) long term correlations may be modeled differently, since 
there are numerous distinct ways to generate power law distri- 
butions (Newman, 2005). For example, as numerically demon- 
strated in Gal and Marom (2013), 1// statistics in neuronal 
firing patterns can be generated by assuming global (cooperative) 
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interactions between ion channels (i.e., not through the volt- 
age). BiophysicaUy, the significance of interactions between ion 
channels is still not clear (Naundorf et al, 2006 and Brief 
Communications arising), but other cellular processes that might 
affect excitability on slower timescales clearly exhibit interac- 
tions (e.g., gene regulation networks Davidson and Levin, 2005). 
Mathematically, such interactions render the slow dynamics 
(Equation 6) non-linear at constant voltage (Gillespie, 2000). It 
would be interesting to generalize the theory we presented here 
in order to understand how to tune the PSD in such a non-linear 
setting, since this has the potential to further reduce the number 
of parameters and model complexity. 

4.2. BIOPHYSICAL IMPLEMENTATION 

We examine our theoretical predictions numerically. We do this 
using a stochastic Hodgkin Huxley type model with slow sodium 
inactivation that was previously fitted to the basic experimen- 
tal results (Soudry and Meir, 2012). We extend this model to 
include four additional slow processes, which resemble slow 
sodium inactivation (Appendix C.2). The only difference is 
that each process is slower than the previous one, and has a 
lower number of ion channels, according to the specific scal- 
ing relation that was derived. The resulting model indeed gen- 
erates f^" noise, and is demonstrated numerically (Figure 3) 
to fit the experimental results of Gal et al. (2010). This is 
the first time, to our knowledge, that a cortical neuron model 
(either biophysical or phenomenological) reproduces experi- 
mental results over such long timescales. Notably, without the 
analytical results, it would be hard to tune the parameters of 
a biophysical neuron, due to the large number of unknown 
parameters. 

Previous works (Lowen et al., 1999; Soen and Braun, 2000) 
demonstrated numerically that, even with constant current stim- 
ulation, incorporating slow processes into an excitable cell model 
can generate in its response. In Lowen et al. (1999) a 
HH model was extended to include multiple slow processes 
with scaled rates in the potassium channel produced f^" firing 
rate response. Their model produced an exponent of a ~ 0.5, 
replicating experiments measurements from the auditory nerves. 
Another work (Soen and Braun, 2000), aiming to reproduce 
the activity of heart cells, produced long term correlations with 
a ^ \.6 — 2 using a reflected diffusion process. 

The identity of the specific slow processes involved in gener- 
ating remains a mystery at this point, since there are many 
possible mechanisms which can modulate the excitability of the 
cell in such long timescales. For example, ion channel numbers, 
conductances and kinetics are constandy being regulated and 
may change over time (e.g., Levitan, 1994; Staub et al., 1997). 
Also, the ionic concentrations in the cell depend on the activ- 
ity of the ionic pumps, which can be affected by metabolism 
(Silver et al., 1997). Finally, the spike initiation region can sig- 
nificantly shift its location with time (e.g., 17 /xm distally during 
48 h of high activity Grubb and Burrone, 2010), and so can cel- 
lular neurites (Paola et al., 2006; Nishiyama et al, 2007). Only 
after such slow processes are quantitatively characterized, can 
we determine their effect on the neuron's excitability at long 
timescales. 



4.3. THE INPUT-OUTPUT RELATION 

The linearized input-output relation of the fitted CBM was 
derived using the methods described in Soudry and Meir (2014). 
This linearized relation decomposes the contributions of external 
inputs and internal fluctuations to the response of the neuron. 
This decomposition (Equation 33) shows that even though the 
neuron can "remember" its intrinsic fluctuations over timescales 
of days, its memory of past pulse inputs can be limited to a shorter 
timescale of ~ 10^ s (Figure 4). Notably, the neuron has this lim- 
ited memory for such inputs even though processes on much 
slower timescales exist in the model. 

In the introduction we mentioned previous works (Lundstrom 
et al., 2008, 2010; Pozzorini et al, 2013) which also described 
the temporal integration in the neuron using power-law filters, 
although in a rather different (non-sparse) stimulation regime. 
Our fitted model indicates that similar power-law integration still 
occurs at very long timescales. However, it is not the input that 
is being integrated, but the internal fluctuations in the neuron, 
and this is what drives the f^" statistics measured by Gal et al. 
(2010). Also, as in Lundstrom et al. (2008, 2010) and Pozzorini 
et al. (2013), the neuronal response in our model is indeed 
affected by the last 10 s of its external inputs. However, our model 
suggests the response will not be significantly affected by spike 
perturbations in its input that occurred more than 10^ s ago. 

Qualitatively, this specific timescale of the input memory stems 
from the "fastest slow negative feedback process" in the model 
(in this specific model, slow sodium inactivation). This process 
responds to perturbations in the input which change the fir- 
ing rate much more quickly than all the other slow processes. 
Its response to perturbation brings the firing rate back to its 
steady state, before slower processes even register that the firing 
rate has changed. Therefore, effectively, these slower processes 
do not store much information about input perturbations. We 
suggest experiments to test input memory directly, by using f^" 
stimulation (Figure 5A), "sum of sines" stimulation (Figure 5B) 
and a variation of the mean stimulation rate (Equation 44 and 
Supplementary Figure 2). 

4.4. SIGNIFICANCE 

This work makes several practical contributions. First, our results 
impose specific constraints on the slow processes that modu- 
late the excitability on very long timescales (e.g., a ratio between 
timescales and channel numbers). Such constraints facilitate the 
construction of neuronal models with "realistic" input-output 
relations over extended timescales. Hopefully, these constraints 
may also help to identify the relevant slow biophysical processes. 
Second, our results suggest that for sparse spiking inputs, the 
memory of a cortical neuron stretches back to the last minute 
of its input, but not much more. This limit could be especially 
relevant when fitting statistical models to neuronal data, and for 
setting limitations on neuronal computations. 

As for the functional significance, it is still not clear why 
the neuronal response fluctuates so wildly, especially at long 
timescales. We end this paper by offering some speculations on 
this issue. We see three possible scenarios. One possibility is 
that these fluctuations are beneficial. For example, such non- 
stationary fluctuations should increase network heterogeneity. 
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which may be advantageous (Tessone et al, 2006; Padmanabhan 
and Urban, 2010). Another possibihty is that these fluctuations do 
not affect neural response when the neuron is connected within a 
network. For example, this could be due to network feedback can- 
celing slow changes in excitability. Finally, it is possible that such 
slow fluctuations are deleterious, an unavoidable "noise" gener- 
ated by the non-stationary environment of the cell. Interestingly, 
f^" noise imposes important constraints on electronic circuits, 
and was predicted to impose similar constraints on neural circuits 
(Sarpeshkar, 1998). 
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